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Abstract 

A general purpose orthotropic elasto-plastic computational constitutive material model has been 
developed to improve predictions of the response of composites subjected to high velocity 
impact. The three-dimensional orthotropic elasto-plastic composite material model is being 
implemented initially for solid elements in LS-DYNA as MATH 3. In order to accurately 
represent the response of a composite, experimental stress-strain curves are utilized as input, 
allowing for a more general material model that can be used on a variety of composite 
applications. The theoretical details are discussed in a companion paper. This paper documents 
the implementation, verification and qualitative validation of the material model using the T800- 
F3900 fiber/resin composite material. 


Introduction 

Composite materials are now beginning to provide uses hitherto reserved for metals in structural 
systems such as airframes and engine containment systems, wraps for repair and rehabilitation 
and ballistic/blast mitigation systems. While material models exist that can be used to simulate 
the response of the materials in these demanding structural systems, they are often designed for 
use with specific materials such as metals [1] [2], polymers [3] and wood [4], Most are tailored 
specifically for an application and have limitations - purely elastic, no rate sensitivity, 
implementation for solid elements only, and limited damage and failure characterization [5]-[10]. 

A few examples of different approaches to composite damage modeling include Littell et al. 
[11], Matzenmiller et al. [12], and others [13]-[18]. Yen’s comprehensive model [19] is 
implemented in LS-DYNA as *MAT_162. Availability of constitutive models for composite 
materials in LS-DYNA is discussed in some length in our companion paper [20] and that 
explains the motivation for this study. 

In the next section we briefly discuss the theoretical development of the constitutive model. This 
is followed by details of the numerical algorithm, verification tests for a composite used in the 
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aerospace industry, the use of the composite in a high speed contact example and finally, 
concluding remarks including details of our ongoing and future work. 

Orthotropic 3D Elasto-Plastic Composite Material Model 


The material law that the model is built upon describes the elastic and permanent deformation of 
the composite with full three-dimensional implementation suitable for solid elements. Only the 
relevant equations used in the numerical algorithm are presented here. The theoretical basis and 
discussions are available in a companion paper [20]. Current development of the model includes 
complete elasto-plastic deformation model, with damage and failure to be added later. 


The Tsai-Wu surface is used as the yield surface for the plasticity model and is defined as 
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The stress components of the yield function coefficients correspond to the yield stresses 
associated with each experimental curve, which vary with the effective plastic strain, thus 
allowing the model to describe different hardening properties in each direction. The non- 
associative flow surface is defined as 

h — 22^22 ^ 33^33 ^^ 12 ^ 11^22 ^^ 23 ^ 22^33 ^^ 31 ^ 33^11 ^44^12 ^55^2i ^ 66^31 


where the H.j terms are the constant flow rule coefficients. The rest of the details can be found 
in the companion paper [20]. 


The first six flow rule coefficients are computed directly from the assumed flow rule coefficient 
value, //jj and the plastic Poisson’s ratios. 
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The last three flow rule coeffieients {Haa, H55, //ee) are calculated iteratively with the user- 
defined range to best fit the shear curves [20]. For example, to determine Haa, the following 
optimization problem is posed using the input curve for the “ 1 - 2 ” shear case. 


Find 7/44 to minimize 


/(ff 44 ) = Z (^22). -(^12). 


(3b) 


such that 

HZ<H,,<HZ (3c) 

where n is the number of data points on the master curve, [ 0 - 22 ), is the effective stress value 

from the master curve, and (<^ 12 ). is the effective stress value for the shear curve using an 
assumed value of Haa. 


Numerical Algorithm 


In this initial implementation, we do not account for rate and temperature dependence, damage 
accumulation nor failure. The following experimental data are needed as LS-DYNA input. 

(1) A total of 12 true stress versus true strain curves under quasi-static test conditions from 
(a) uniaxial tension in 1, 2 and 3-directions, (b) uniaxial compression in 1, 2 and 3- 
directions, (c) shear in 1-2, 2-3 and 3-1 planes and (d) off-axis (e.g. 45 degrees) uniaxial 
tension or compression in 1-2, 2-3 and 3-1 planes. The tabulated x-y data for each curve 
is read in and used appropriately. 

(2) The modulus of elasticity, Poisson’s ratio and average plastic Poisson’s ratio obtained 
from the tension and compression tests. 


The effective plastic strain can be written in terms of the experimental stress versus total strain 
data. In the 1 -direction for example, this is written as 






<i(<) 


( 4 ) 


where <t[j is the experimental compressive true stress in the 1 -direction, is the total true 
strain in the 1 -direction, £'jj is the elastic modulus in the 1 -direction, is the true plastic strain 
in the x-direction, is the effective plastic strain and h is the flow rule value, which is 
equivalent to the effective stress, shown in Eqn. (2). The material model uses a typical elastic 
stress update as 

=o"+CAt:(£-£^) (5) 

where C is the ortho tropic stiffness matrix. At is the time step, £ is the total strain rate and £^ 
is the plastic strain rate. The stiffness matrix is written in terms of the compliance matrix as 
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where are the elastie moduli in the principal material directions, G-j are the elastic shear 
moduli and v-j are the elastic Poisson’s ratios. The yield stresses used to determine the flow rule 
coefficients are summarized into a single vector, corresponding to each of the experimental test 
curves as q ctjj ct 33 ctjj ctjj ct 33 cTjj ctjj ct 3 j <^4 s-u *^ 45-23 *^ 45 - 3 i] • Lastly, the 

consistency condition is written in terms of the rate of yield function as 

(V) 
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which is defined as a non-associative flow with strain hardening (hence the inclusion of the yield 
stress vector). Eqn. (7) can be expanded as 
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where d is written in terms of the constitutive matrix and strain rates and X is the effective 
plastic strain. Solving for the effective plastic strain rate produces the following equation. 
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Initially, a perfectly elastic response is assumed. Thus, an elastic predictor is used to compute the 
stresses as 

(cTii) =(cTjj) +CjjAt(^jj) + Cj2At(^22) + L'[3At(£33) 
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With the elastic trial stresses computed, a trial yield function value can be calculated and used to 
determine if the load step is elastic or plastic. It should be noted that the elastic trial stress is used 
as the stress in the next step if the trial yield function value is less than or equal to zero. 
Otherwise, the stress state is beyond the yield surface and must be brought back using a radial 
return and plastic corrector, AA . If the trial yield function happens to be greater than zero, then 
AA must be greater than zero. The value of AA is determined using a secant iteration, with 
A/l' = 0 for the first iteration and the second value for the increment of the plastic multiplier 
determined from consistency, 

df 
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where the derivatives of q are taken as zero, meaning that the response is perfectly plastic. The 
initial calculation of the second increment corresponds to perfect plasticity in order to ensure that 
the stress state returns to the interior of the yield surface, resulting in a negative value of the 
yield function and bounding the solution. From each estimate of AA , the corresponding 
corrected stresses can be computed as well as the yield stresses obtained from the input stress 
versus effective plastic strain (A) data. With these terms calculated for both estimates of AA , 
the yield function value for each can be computed with the third estimate evaluated as 

AA^=AA^-f^j^^ (12) 

The corrected stresses and yield stresses associated with AA^ are then calculated and used to 
obtain a new estimate of the yield function value. Convergence of the secant iteration is 
determined by the following conditions 

/3«o^a^=a^^ (13) 


AA" = AA" 


AA" = AA" 


with the secant iteration continuing if the new yield function value is not less than some defined 
tolerance. Once convergence is met, the plastic multiplier increment to return the stress state to 
the yield surface is known and the stresses can be updated as 


_«+l 
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(14) 


Finally, the yield stresses are updated as well, using the new value of the overall effective plastic 
strain, A , in each input curve to determine the corresponding yield stress level as 

=^(r +A;C) (15) 

This results in anisotropic strain hardening, as a yield stress increase in each direction is initiated 
with plasticity in any direction, but at different levels. A detailed algorithm is presented below. 


Step 1: Preprocessing 

(a) Convert stress/strain input curves to stress versus effective plastic strain using Eqn. (4). 

(b) Store initial yield stresses in q . 

(c) Calculate and store the elastic moduli from the initial yield stress and strain values. 

(d) Store the three elastic Poisson’s ratio values and flow rule coefficients from input. 

(e) Compute optimal values of the flow rule coefficients so as to match the input curves as 
closely as possible. 
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Step 2: Initialization 

Parameters available at start of step: q„ , \ . 


Step 3: Elastie predietor 

Compute and set the yield funetion eoeffieients using Eqn. (la). 

Construct the constitutive matrix using Eqn. (7). 

Compute elastic trial stresses, , using Eqn. (1 1). 

Compute the trial yield function, , using the elastic trial stresses in Eqn. (1). 

If <0 , the elastic solution is correct, set A/l,^ = 0 and go to stress update. Else go to 
plastic corrector. 


Step 4: Plastic corrector 
Set A/l' = 0 . 


Calculate AA^ from Eqn. (11). 

Eoop through secant iteration for specified max number of iterations: 

Compute the new estimate of the stress for each plastic multiplier increment ^AA*, A/l^^ 
using Eqn. (14). 

Calculate the effective plastic strains at the next time step using the estimates of 


the plastic multiplier increment; then obtain the strain hardening parameters with the 
effective plastic strain from the curve data (true plastic stress as a function of effective 
plastic strain). 


Compute new plastic multiplier increment, . Calculate new estimate for the yield 
function, .If f^< tol , set A/l = A/l^ , exit the loop and go to stress update. Else update 


secant iteration parameters using Eqn. (14) and go to next step of secant iteration. 


Step 5: Stress Update 

Calculate using Eqn. (14). 


Step 6: Update history variables for plastic work and work hardening parameters ( q and A ) 

Set 4+1 =4 +M,- 

Determine new yield stresses, q„^j, using 4+i ^^d Eqn. (15). 


Model Verification 

The composite material model was tested and verified using experimental data obtained from 
T800S/3900-2B [P2352W-19] BMS8-276 Rev-H-Unitape fiber/resin unidirectional composite 
[21]. Toray describes T800S as an intermediate modulus, high tensile strength graphite fiber. The 
epoxy resin system is labeled E3900 where a toughened epoxy is combined with small 
elastomeric particles to form a compliant interface or interleaf between fiber plies to resist 
impact damage and delamination [22]. Magnified views of the composite are shown in Eigs. 1 
and 2. 
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Fig. 1. Side view obtained using optical microscopy 



Fig. 2. Logitudinal view obtained using a SEM 


The source of the input data for the model (both experimental and virtual) is listed Table 1. In the 
Data column, Experimental refers to experimental data generated at Wichita State [21], and 
MAC-GMC [23] and VTSS [24] refer to the use of numerical simulation techniques to generate 
the stress-strain curve when experimental data are not available. MAC-GMC uses the generalized 
method of cells whereas VTSS uses displacement-based finite element analysis in the context of 
virtual testing to predict elastic and inelastic material properties. 


Table 1. Generation of input data for T800-F3900 composite 


Curve 

Data 

Tension Test (1 -Direction) 

Experimental 

Tension Test (2-Direction) 

MAC-GMC and VTSS 

Tension Test (3 -Direction) 

Transverse isotropy 

Compression Test (1 -Direction) 

Experimental 

Compression Test (2-Direction) 

MAC-GMC and VTSS 

Compression Test (3 -Direction) 

Transverse isotropy 

Pure Shear Test (1-2 Plane) 

Experimental 

Pure Shear Test (2-3 Plane) 

MAC-GMC and VTSS 

Pure Shear Test (1-3 Plane) 

Transverse isotropy 

Off-Axis Test (45°, 1-2 Plane) 

MAC-GMC and VTSS 

Off-Axis Test (45°, 2-3 Plane) 

MAC-GMC and VTSS 

Off-Axis Test (45°, 1-3 Plane) 

Transverse isotropy 


The fiber (transversely isotropic, linear elastic) and the matrix (isotropic, elastic perfectly plastic) 
properties are listed in Table 2. The properties for the latter are not publicly available. Hence 
they were assumed as shown in the table. 

Table 2. T800 -F3900 fiber and matrix properties (Volume Fr action = 0.54) 


Property 

Fiber 

Matrix 

El, psi 

4(10’) 

5(10^) 

E 2 , psi 

2.25(10’) 


E3, psi 

2.25(10’) 


V12 

0.2 

0.35 

V23 

0.25 


V13 

0.25 


G, psi 

1.5(10’) 

1.85(10^) 

ay, psi 


2(10^) 
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Results from the preprocessing step in terms of computing the optimal flow-rule coefficients 
(solution to Eqns. (3b) and (3c)) are shown in Figs. 3 and 4. 


Optimized H44 Value-Generated Curve 



Fig. 3. Comparison of experimental curve with optimized H 44 (and Hee) value curve 


Optimized Hss Value-Generated Curve 



Fig. 4. Comparison of experimental curve with optimized Hss value curve 

With the computed data so far, the verification process can be executed. In the verification test 
results shown in the figures that follow, the term Experimental (Input) refers to experimental 
data (when available) or to data created via virtual testing, and the term Simulated (x-element) 
refers to the data created using results from finite element simulations where x is the number of 
finite elements in the FE model. 

Schematics for the tension and compression tests are shown in Fig. 5 and Fig. 6, respectively. 
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Fig. 5. Schematics for tension test simulation cases (a) 1 -direction (b) 2 and 3-directions 


2(y) 


1 (y) 
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(b) 


Fig. 6. Schematics for compression test simulation cases (a) 1-direction (b) 2 and 3-directions 

The simulations for the tension and eompression tests in all three direetions were exeeuted using 
different mesh sizes ranging from 1 element to 64 elements. 


64-element meshes for the tension and eompression test eases are shown in Fig. 7 and Fig. 8, 
respeetively. 



Fig. 7. 64-element mesh for tension test cases 


Fig. 8. 64-element mesh for compression test cases 


For both tension and eompression tests, nodes on faee ABC were eonstrained in the x-direetion 
and the eenter node was also eonstrained in the y-direetion. A 0.5 in/s displaeement in the x- 
direetion was applied to all the nodes on the right faee, E-D. For the eompression tests, eare was 
taken to ensure that the model yielded but did not buekle. The simulated and experimental stress- 
strain eurves for the tension and eompression tests are shown in Figs. 9 and 10. 
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Fig. 9. Simulated and experimental stress-strain curves for 1-direction (a) tension and (b) compression 


90-Degree Tension Test (2/3-direction) 
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90-Degree Compression Test (2/3-direction) 



(b) 

Fig. 10. Simulated and experimental stress-strain curves for 2/3-directions (a) tension and (b) compression 


A schematic for the pure shear test is shown in Fig. 11. Nodes on the bottom surface are 
restrained in the x and y-directions while nodes on the top surface are restrained in the y- 
direction. A displacement controlled loading of 0.5 in/s in the x-direction is applied to all the 
nodes on the top surface. 



V 


3(z) 

Fig. 11. Schematic for pure shear test in 1-2/3-1 Fig. 12. 64-element mesh for pure shear test cases 
plane 

The 64-element mesh is shown in Fig. 12. The simulated and experimental stress-strain curves in 
the 1-2/3-1 plane are shown in Fig. 13 and in the 2-3 plane in Fig. 14. 
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1-2/3-1 Shear Test 



Fig. 13. Simulated and experimental stress-strain curves for pure shear in the 1 -2/3-1 plane 

2-3 Shear Test 



Fig. 14. Simulated and experimental stress-strain curves for pure shear in the 2-3 plane 


The (tension) off-axis test model is shown in Fig. 15. Nodes on face ABC were constrained in 
the x-direction and the center node was also constrained in the y-direction. The simulated and 
experimental stress-strain curves are compared in Fig. 16, showing very little differences in the 
results. Finally, Fig. 17 shows the results from a loading-unloading-reloading test case. 
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Fig. 15. Schematic for 45° off axis test in 1-2/3-1 plane 
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Fig. 16. Simulated and experimental stress-strain curves for 45° off axis test in 1-2/3-1 plane 


Unloading/Reloading 90 Degree Test (64-Element model) 
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Fig. 17. Simulated and experimental stress-strain curves for unloading/reloading in the 2-direction 
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Simulations of all the 12 curves as an integral part of the verification process shows very good 
correlation with the input data. We have used boundary conditions that resemble the 
experimental test setup as closely as possible so as to capture the actual test conditions. The use 
of 4 different meshes with increasing mesh density, shows that the results are mesh independent. 


Structural Test 

To gage the performance of the implemented constitutive model, a structural test case is 
designed to qualitatively validate the developed model. A simply-supported 12” long x 2.4” wide 
X 0.6” thick plate made of T800-F3900 composite material is subjected to an impact from a steel 
ball. Fig. 18 shows the schematic diagram. The nodes on the bottom of the left edge are 
constrained in the x, y and z-directions while nodes on the bottom of the right edge were 
constrained only in the y-direction. In the absence of experimental data, the primary objectives 
were to ensure that the material model (in the absence of damage and failure) would yield 
expected results. That is, the ball would impact the plate and bounce back, and the peak stresses 
in the plate would occur around the point of contact immediately after the impact. The LS- 
DYNA FE model is shown in Fig. 19. References to LS-DYNA keywords can be found in 
Hallquist [25]. 

-^ 1 . 20 " 


y 

b 

o 



Fig. 18. Schematic diagram of the test 




Fig. 19. LS-DYNA FE model at the start of the 
analysis 



The model has 6,478 nodes and 5,707 elements of which 1,147 8-noded hexahedral elements are 
used to model the plate (master surface) and the remaining 4,560 elements are used in the ball 
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(slave surface). The ball is modeled using *MAT_3 and the plate as *MAT_213. 
*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE is used to describe the contact 
between the ball and plate. Both the static (fs) and the dynamic (fd) coefficients of friction are set 
to 0.1. A 2% viscous damping coefficient is used. Segment-based contact is used by setting the 
soft parameter to 2. The steel ball is given a -700 in/s initial velocity in the y-direction. The 
analysis is carried out for 0.1 seconds. The results are as expected. The exaggerated deformed 
shape is shown in Eig. 20. The von Mises stress throughout the test duration is plotted for the 
elements closest to the point of contact in Eig. 21. 



Time (s) 


Fig. 21. von Mises stress vs. time for elements near point of impact 


Conclusions 

A generalized three-dimensional composite material model is developed and verified. The 
preliminary results are encouraging as the implemented constitutive model is able to reproduce 
the experimental stress-strain curves. In addition, a simple contact problem involving the 
composite material yields reasonable results. Euture work includes the following: adding damage 
accumulation and failure in the material model, supporting temperature and strain rate dependent 
stress-strain curves, using alternate methods to compute the plastic multiplier when the initial 
yield and flow surface convexity conditions are not met, using a fiber-referenced convected 
system to track the principal material directions and supporting plate/shell elements and to 
validate the developed model using data from several structural tests. 
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